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(57) Abstract 

A process is disclosed for predicting a characteristic of a fluid-containing physical system such as a reservoir. The physical system 
is discretized to a plurality of gridcells. For each gridcell, nonlinear governing equations are constructed that are representative of fluid 
properties for each fluid. The nonlinear terms are linearized to derive linear equations at a current timestep. The linear equations are then 
solved to determine estimates of the solution. The estimated solution is then improved by a gridcell-by-gridcell computation wherein flows 
into and out of a particular gridcell are determined using linearized values of nonlinear terms in gridcells adjacent that particular gridcell. 
The improved solution is used to predict a property of a fluid in the physical system at end of the timestep. These calculation steps are 
repeated for a plurality of timesteps and the results are used to predict a property of the physical system and the fluids it contains as a 
function of time. 
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IMPROVED PROCESS FOR PREDICTING BEHAVIOR OF A 
SUBTERRANEAN FORMATION 



5 This application claims the benefit of United States Provisional Application 

No. 60/074188 filed February 10, 1998. 

FIELD OF THE INVENTION 

The invention relates generally to exploitation and development of hydrocarbons 
in an underground reservoir and, more particularly, to an improved process for 
1 0 predicting the behavior of a subterranean, hydrocarbon-bearing formation. 

BACKGROUND OF THE INVENTION 

Reservoir simulation is a process of inferring the behavior of a real reservoir 
from the performance of a model of that reservoir. Because mass transfer and fluTd flow 
15 processes in petroleum reservoirs are so complex, reservoir simulations can only be done 
using computers. Computer programs that perform calculations to simulate reservoirs 
are called reservoir simulators. The objective of reservoir simulation is to understand the 
complex chemical, physical, and fluid flow processes occurring in a petroleum reservoir 
sufficiently well to be able to predict future behavior of a reservoir and to maximize 
20 recovery of hydrocarbons. The reservoir simulator can solve reservoir problems that are 
not solvable in any other way. For example, a reservoir simulator can predict the 
consequences of reservoir management decisions. Reservoir simulation often refers to 
the hydrodynamics of flow within a reservoir, but in a larger sense it also refers to the 
total petroleum system which includes the reservoir, the surface facilities, and any 
25 interrelated significant activity. 

Figure 1 illustrates schematically four basic steps in conventional reservoir 
simulation of a petroleum reservoir. The first step (step 1 in Fig. 1) is to construct a 
mathematical model of a real reservoir based on the chemical, physical, and fluid flow 
processes occurring in the reservoir. The mathematical model will comprise a set of 
30 nonlinear partial differential equations. The second step (step 2 in Fig. 1) involves 
discretization of the reservoir in both time and space. Space is discretized by dividing 
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the reservoir into suitable cells with each cell having a set of nonlinear finite difference 
equations. The third step (step 3 in Fig. 1) is to linearize the nonlinear terms that appear 
in the nonlinear finite difference equations and, based on this linearization, construct 
linear algebraic equations. The fourth step (step 4 in Fig. 1) is to solve the linear 
5 algebraic equations. The solution provides a prediction of reservoir behavior, which 
enables a petroleum engineer to predict reservoir performance, including the rate at 
which the reservoir can be produced. The accuracy of the model can be checked against 
the history of the reservoir after the model has been subjected to a simulated recovery 
process. 

10 In developing a mathematical model (step 1 in Fig. 1), the following principles 

apply: 

(1) Anything that enters or leaves the reservoir must cross a boundary. 

(2) At some initial time the reservoir can be described by some set of conditions. 

(3) The processes that occur within the reservoir obey known physical laws and 
1 5 consequently can be described by a set of conditions. 

From these principles, a complete mathematical model of the reservoir can be developed 
which is a combination of (i) boundary conditions, (ii) initial conditions to describe 
conditions at zero time, and (iii) a system of equations that govern the behavior of the 
fluids in the reservoir. The governing equations are obtained by combining the 
20 following physical principles: 

(1) conservation of mass, 

(2) conservation of momentum, 

(3) conservation of energy (first law of thermodynamics), and 

(4) thermodynamic equilibrium. 

25 The resulting equations are nonlinear partial differential equations that describe the 
unsteady-state flow of all fluids throughout the reservoir and relate the pressure and 
saturation changes of the fluids with time throughout the reservoir. Examples of 
methods for constructing mathematical models of reservoirs are described in the 
following books: Peaceman, D. W., Fundamentals of Numerical Reservoir Simulation, 

30 Elsevier Scientific Publishing Company, Amsterdam, (1977); and Mattax, C. C. and 
Dalton, R. L., Reservoir Simulation^ Monograph Volume 13, Society of Petroleum 
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Engineers (1990). 

Solution of the mathematical model requires determining values of the 
independent parameters that satisfy all the governing equations and boundary conditions. 
In general, there are two solution choices: analytical or numerical methods. Analytical 
5 methods are not used because the equations are highly nonlinear and are practically 
impossible to solve by today's methods. The most common method for solving the 
equations uses a discretization process. 

The second step of conventional reservoir simulation (step 2 in Fig. 1) involves 
subdivision of distance and time into definite, specified increments. This discretization 

10 process transforms the continuous partial differential equations to a finite dimensional 
system of algebraic equations (commonly called finite difference equations). Space is 
discretized by subdividing the reservoir volume into contiguous cells conforming to a 
predetermined grid pattern, frequently called "gridblocks" or "gridcells." A typical 
reservoir model can have thousands of cells. These cells will have sides with lengths 

15 ranging from a few meters to a few hundred meters. 

The life of a reservoir is discretized (or divided) into time increments. A 
simulator computes changes in each gridcell (flow, pressure, etc.) over each of many 
timesteps. Typically, conditions are defined only at the beginning and end of a timestep, 
and nothing is defined at any intermediate time within a timestep. Consequently, 

20 conditions within each gridcell may change abruptly from one timestep to the next. 

Usually, timesteps are chosen to be small enough to limit sizes of these abrupt changes 
to acceptable limits. The size of the timesteps depends on accuracy considerations and 
stability constraints; generally the smaller the timestep the more accurate the solution, 
but smaller timesteps require more computational work. 

25 Although discretization of the reservoir is an abstraction, it is qualitatively 

correct to visualize gridcells as well-stirred tanks with permeable sides. The contents of 
a gridcell are considered uniformly distributed within the cell and the rates at which 
fluids flow in or out are determined by the permeabilities of the sides of the cell and the 
pressure differences between adjacent gridcells. In essence, the mathematical problem is 

30 reduced to a calculation of flow between adjacent gridcells. 

Because of the discretization process, discontinuities in saturation and pressure 
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exist at interfaces between gridcells. Similarly, saturation and pressure are also 
discontinuous in time, with the reservoir condition defined only at the end of each 
timestep. Properties that are functions of saturation (such as relative permeability and 
capillary pressure) can change significantly over a timestep. In the real reservoir, at any 
5 particular time, a gridcell has only a single value of each phase saturation and any 

property that is saturation dependent. To represent variations in reservoir properties, the 
properties of the gridcells must differ from cell to cell. The abruptness with which a 
property changes between neighboring cells is largely a function of gridcell size. During 
a timestep, the saturation-dependent properties will also have assumed many values. 

10 However, in reservoir simulation a single value for each saturation-dependent variable 
must be used to describe conditions over the entire timestep. It is therefore apparent that 
the finite difference equations can only approximate the partial differential equations that 
hold for any point within the grid system at a predetermined timestep. 

The third step in conventional reservoir simulation (step 3 in Fig. 1) is to 

1 5 linearize the nonlinear terms that appear in the nonlinear finite difference equations and, 
based on this linearization, construct a linear set of algebraic equations. A linearization 
represents a term as a straight line function of the quantity or quantities upon which it 
depends. 

The fourth step in conventional reservoir simulation (step 4 in Fig. 1) is to solve 
20 the algebraic equations. It is necessary to solve the equations for the unknown values, 
which in a reservoir containing oil and water will normally include two of the following 
four variables: oil pressure, water pressure, oil saturation, and water saturation. Other 
quantities can be derived from these variables, such as oil production rate and water 
production rate. 

25 Many simulation methods have been proposed. The method chosen can affect 

the stability and accuracy of the solution and some methods require more computational 
work than other methods on a per-timestep basis. The methods differ primarily on how 
they treat the way the reservoir variables (such as pressure and saturation) vary in time. 
Most methods involve variations of the following two procedures: 

30 (1) Explicit procedures use mobilities and capillary pressures computed as functions 
of saturations at the beginning of a timestep. The saturations are known from the 
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previous timestep calculations. The mobilities and capillary pressures are 
assumed to maintain the same values during a timestep that they had at the 
beginning of the timestep. 
(2) Implicit procedures use mobility and capillary pressure calculated as functions of 
5 saturation at the end of the timestep. The values are not known until calculations 

for the timestep have been completed. As a result, they must be determined 
using an iterative process. 

A brief overview of six basic methods for solving finite difference equations is 

given below. These are the Fully Explicit method, the Implicit Pressure Explicit 
1 0 Saturation method (IMPES), the Fully Implicit method, the Sequential Implicit method 

(SEQ), the Adaptive Implicit method (AIM), and the Cascade method. 

Fully Explicit Method 

The Fully Explicit method is the simplest solution method because for the 

purposes of computing flow between cells it assumes that pressures and saturations 
1 5 remain constant at their initial values throughout each timestep. In effect, the inter-cell 

flows for each timestep are computed without taking into account the way pressures and 

saturation are changing. As a result, this method is prone to instability. For example, 

the procedure may predict that the total mass of a fluid within a cell becomes negative. 

As the timestep increases, the probability of instability increases until, at some 
20 sufficiently long timestep, unphysical predictions are inevitable. In reservoir simulation, 

these considerations limit the timestep to such a small value that the Fully Explicit 

method is of no practical use. 

IMPES Method 

In the IMPES method, saturations are assumed to remain constant at their initial 
25 values throughout the timestep. Pressures, on the other hand, are assumed to reach their 
final values very quickly and to remain there throughout the timestep. In this method, 
the flow rates and the pressures at the end of the timestep are interdependent and must be 
determined by iteration. The term "implicit*' is used because of the interdependence of 
the flow rate during a timestep and the pressure solution at the end of the timestep. The 
30 pressure- and saturation-dependent terms lag behind the pressure terms. The basic 
procedure is to obtain a single pressure equation by a combination of flow equations. 
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After the pressure has been advanced in time, the saturations are updated explicitly. 
When the saturations are calculated, new capillary pressures are calculated, which are 
explicitly used at the next timestep. Implied in this approach is the belief (or hope) that 
the saturation and pressure values are not changing very rapidly. In some models the 
5 method produces good results. However, the method becomes unstable when there is 
rapid fluid movement in the reservoir (around the wells, for example) or when capillary 
pressures affect flow strongly. For this reason, IMPES is often not practical. 
Fully Implicit Method 

The Fully Implicit method, which is more complex than the Fully Explicit and 

10 IMPES methods, is unconditionally stable because it treats both pressure and saturations 
implicitly. Flow rates are computed using phase pressures and saturations at the end of 
each timestep. In this method, saturations cannot fall below zero because a fluid can 
flow only if it is mobile at the at the end of a timestep, and fluids are mobile only for 
saturations greater than zero. Saturations cannot exceed unity because this would imply 

15 negative saturations in other phases. The calculation of flow rates and pressure and 

saturation solutions involves the solution of nonlinear equations using a suitable iterative 
technique. As the pressures and saturations are solved, the updating of these terms 
continues using new values of pressure and saturation. The iteration process terminates 
when the convergence criteria are satisfied. 

20 The Fully Implicit method is stable for all timestep lengths, and guarantees that 

saturations remain within physical bounds. The unconditional stability of the Fully 
Implicit method makes it possible for the simulation to take much longer timesteps than 
would be possible for a non-implicit computation. Although the Fully Implicit method 
requires an order of magnitude more computing time per timestep than the IMPES 

25 method, the timestep can be several orders of magnitude longer, resulting in a net 
decrease in computing time. 

The main drawback of the Fully Implicit method is the amount of computer time 
that it requires. In terms of computing cost, the method is generally satisfactory in 
models of single wells or parts of a reservoir, but it can be quite expensive to use in 

30 models of entire reservoirs. Several attempts have been made to reduce the 

computations required, possibly at the cost of accepting a method that does not permit 
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the timestep sizes of the Fully Implicit method. The sequential implicit method, the 
adaptive implicit method, and the Cascade method have been proposed as ways of 
reducing the computational time. 
Sequential Implicit Method 
5 The Sequential Implicit method (SEQ method) incorporates implicit treatment of 

saturations, but without solving simultaneously for pressures and saturations. The SEQ 
method consists of two steps. The first step solves a set of pressure equations in exactly 
the same way as is done in the IMPES method. This set comprises a single equation per 
cell, and solving it yields a complete, new pressure distribution at the end of a timestep. 

10 In a second step, the pressure distribution is used to calculate the sum of the velocities of 
all phases at each boundary between cells. These total velocities are used in constructing 
a set of saturation equations. This set comprises two equations per cell in three phase 
cases and one equation per cell in two phase cases and is solved simultaneously to yield 
saturations at the new time. The second step is an implicit solution for saturations using 

15 linearized implicit velocities. Saturations in each cell are determined by using implicit 
(end-of-time-step) values of relative permeabilities in inter-cell fluid flow terms. This 
method requires simultaneous solution of all saturation equations. The per-timestep 
computational cost of the SEQ method, though 50 to 100 percent greater than that of the 
IMPES method, is typically about one-fifth that of the Fully Implicit method. The SEQ 

20 method can take much longer timesteps than the IMPES method, and in some models its 
timestep size approaches that of the Fully Implicit method; in other models, the SEQ 
method requires up to several times as many timesteps as the Fully Implicit method. 
Therefore, in some models the SEQ method takes less computing time overall than the 
Fully Implicit method; in other models, it takes more. 

25 Adaptive Implicit Method 

The Adaptive Implicit method (AIM) is sometimes used in an attempt to combine 
the superior timestep size of the Fully Implicit method with the low computational cost 
of IMPES. It is based on the observation that only small regions in the typical model 
require the stability of Fully Implicit computations, while a simpler method would be 

30 adequate in the rest of the model. AIM performs Fully Implicit computations only in the 
parts of the model where they are believed to be needed and uses IMPES in the rest of 
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the model. Typically, Fully Implicit computations are made in only a small part of the 
model, containing perhaps ten percent of the model's cells. 

The discussion above assumes that the Fully Implicit method can take timesteps 
of arbitrary length. In theory, this is true. Practically, there are difficulties arising from 
5 the need to solve a nonlinear system of equations. The equations are solved using an 
iterative technique called "the Newton method," which makes successive estimates of 
the solution. This procedure is often referred to as a Newton iteration. When the 
Newton method works well, each new estimate is closer to the correct solution than its 
predecessors. Sometimes, it does not work well, and convergence to the correct solution 

10 is obtained slowly or not at all. When convergence is poor, the convergence problem is 
typically confined to a very few cells. Despite this, the full Newton iteration 
computation is normally carried out at all cells (or all implicit cells if AIM is in use). A 
possible way to improve computational efficiency is to, in some manner, focus the 
computational effort on the few cells having difficulty. 

1 5 Cascade Method 

The Cascade method provides a way of improving the computational efficiency. 
This method computes saturation one phase at a time, one cell at a time. At a given 
point in the computation, the cell at which saturation is currently being computed is 
selected to be one for which all incoming flows of the selected phase have already been 

20 computed. Thus, it is only necessary to compute flows of the phase out of the cell. 

Assuming that upstream weighting is being used, these flows depend on saturation in the 
current cell only. 

The first step in the Cascade method is to compute pressure, just like the EMPES 
method. Once all pressures are known, it is always possible to find a cell for which all 

25 incoming flows of a selected phase have already been computed. Initially, this is likely 
to be a cell into which the fluid phase is being injected. Once saturation in this cell has 
been computed, all flows into at least one of its neighbors will be known. Saturation in 
this cell can be computed, after which all flows into at least one more cell will be known, 
etc. In practice, the computation is performed in effect by sorting the cells such that the 

30 one with the largest potential is first, the one with the next largest potential second, and 
so on down to the one with the smallest potential, which is last. Then the cell-by-cell 
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computations are performed in this order. Since all flow into a given cell must come 
from cells having larger potentials, all flows into the current cell are known at any stage 
of the computation. 

The Cascade method achieves its effectiveness by localizing the nonlinear 
5 computations to a single cell and, typically, a single unknown within that cell. Because 
this problem is nonlinear, it must be solved iteratively. Because the computation 
involves only one unknown, an iterative procedure to solve it can be selected for 
reliability rather than efficiency. Only cells at which saturation is changing significantly 
will require more than one iteration, and only cells at which it is changing rapidly will 

10 require more than two or three iterations. As a result, the computational effort is 
concentrated in the cells in which it is needed most. 

The Cascade method is particularly suited for models in which flow of all phases 
is in the same direction. This condition almost never occurs in practical reservoir 
models, since phase, density differences induce countercurrent flow. To some extent, this 

15 problem can be overcome by performing a separate Cascade method computation for 
each phase. However, the method has other drawbacks. First, in the presence of 
capillary forces, the assumption that flow of a phase out of a cell depends only on that 
phase's saturation in the cell is no longer valid. Second, an implicit assumption of the 
Cascade method is that saturation of a given phase is affected only by flow of that phase. 

20 In the presence of interphase mass transfer, this assumption is not satisfied. Perhaps 
because of these shortcomings, the Cascade method is not widely used. 
Linearizations 

The linearization step (step 3 of Fig. 1) and the solution step (step 4 of Fig. 1) 
are dependent on each other. In the process of linearization, the algebraic equations will 

25 have different forms depending on the solution technique chosen. For example, the 
IMPES method linearizes only the pressure-dependent terms, such as specific volume. 
The specific volume is therefore expressed as a linear function of pressure. The SEQ 
method linearizes the same pressure-dependent terms with respect to pressure, and it also 
linearizes phase fractional flow terms with respect to saturations. The Fully Implicit 

30 method linearizes the pressure-dependent terms with respect to pressure and the 

saturation-dependent terms (which comprise the relative permeabilities and capillary 
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pressures) with respect to saturations. 

The Fully Implicit method uses the linearization differently from the IMPES 
method and the SEQ method. In the IMPES and SEQ methods, solving the linearized 
equations gives a solution at the end of each timestep. In the Fully Implicit method, 
5 which uses the Newton method, solving the linearized equations yields approximate 
solutions. Newton iterations are repeated until the resulting estimates of the solutions 
are considered to be accurate enough based on pre-specified convergence criteria. 

One difficulty with conventional linearization procedures is that under certain 
reservoir conditions the computed saturations can be negative, a physical impossibility. 

10 If subsequent calculations are made based on the impossible saturations, erroneous 
results are produced. These erroneous results can occur, for example, in a gridcell 
having a high fluid throughput during a timestep. This occurs frequently near wells, or 
in gridcells that simulate fractures in a reservoir, where gridcells are small and the fluid 
velocities are large. Hundreds or thousands of gridcell pore volumes can flow through a 

15 gridcell in a single timestep. This high throughput can cause stability problemsTin 
reservoir simulation. 

The most common way to deal with this potential instability is to use the Fully 
Implicit method, using a Newton iteration at each timestep with damping used to prevent 
the iteration from diverging. Because a coupled set of equations is being solved, the 

20 Fully Implicit method can be expensive computationally. The SEQ method has also 
been proposed for high throughput gridcells, but it is typically less stable and usually 
takes smaller timesteps than the Fully Implicit method. 

A need exists for an improved reservoir simulation computation method that is 
more reliable and accurate than the computational methods used in the past and at the 

25 same time is computationally efficient. 



SUMMARY 

The process of this invention predicts properties of fluids in a fluid-containing 
physical system, and is particularly suitable for predicting the behavior of a subterranean 
30 formation containing multiple fluids (such as oil, gas, and water) and for predicting 

behavior of multiple fluids in wellbore tubing, surface flowlines, separators and related 
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surface facilities. This simulation process overcomes the shortcomings encountered in 
the past in using the Fully Implicit method and the SEQ method. Compared to 
conventional simulators that use the Fully Implicit method or the SEQ method, this 
invention enables use of longer timesteps when using the SEQ method, and it can reduce 
5 the number of iterations for convergence when using the Full Implicit method. 

The first step in the process is to equate the physical system being simulated to 
a volumetric grid system comprising a plurality of gridcells. For each gridcell, 
nonlinear governing equations are constructed that are representative of rock and fluid 
properties in each gridcell for each fluid. These equations, which represent the fluid 

10 and transport properties as implicit in time, can be constructed using the finite 

difference method, the finite element method, or the finite volume method. These 
fluid and transport properties are nonlinear functions of the solution variables and are 
determined by performing iterative computations. In each iteration of this iterative 
process, the nonlinear terms are linearized to derive linear equations at a current 

1 5 timestep. The linear equations are then solved to determine an estimated solution at 
the end of the time interval. Then, linear estimates of the nonlinear fluid and transport 
properties are computed using the linearizations employed to construct the linear 
equations, evaluating these linearization at the estimated solution. An improved 
solution is then determined in a gridcell-by-gridcell computation. This computation 

20 requires computing flows into and out of a particular gridcell, which in turn requires 
fluid and transport properties at the end of the time interval in the particular gridcell 
and in the gridcells adjacent to it. In computing these flows, the fluid and transport 
properties in the particular gridcell are determined based on the improved solution, 
while the linear estimates of the fluid and transport properties are used for these 

25 properties in one or more gridcells adjacent to the particular gridcell. The improved 
solution is then used to predict a property of one or more fluids in the physical 
system, such as fluid saturation. These calculation steps are repeated for a plurality of 
timesteps and the results are used to predict a property of the physical system and the 
fluids it contains as a function of time. The iteration steps are preferably repeated to 

30 satisfy a predetermined convergence criterion. Once convergence is reached, the 
resulting solution provides a prediction of fluids in the physical system. 
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If the process of this invention is applied to the SEQ method, after computing the 
pressure and "predicted saturations" in the conventional SEQ method, additional 
computations are made to update the linearized fractional flows. The updated linearized 
fractional flows are more accurate than the saturations from which they are computed. 
5 The quantities upon which the updated linearized fractional flows are based, which are 
the updated linearized relative permeabilities and capillary pressures, are then used as 
input to a gridcell-by-gridcell calculation of the saturations. In this calculation, adjacent 
gridcell pressures and saturation-dependent terms are treated as "known," with the 
known values being those resulting from the updated linearization. Stated another way, 
10 a gridcelPs saturations are computed while holding fixed the adjacent gridcells' 
linearized relative permeabilities, capillary pressures, and fractional flows. The 
computed saturations from the additional computation are reasonable and can be used in 
the next outer iteration in the SEQ method. 

If the process of this invention is applied to the Fully Implicit method, the 
15 customary Newton method calculation is used to obtain new estimates of pressure and 
saturation. Linearized relative permeabilities and capillary pressures are then updated 
based on these new estimates of pressure and saturations. The linearized quantities are 
then used as input to a gridcell-by-gridcell calculation of the saturations, wherein 
adjacent gridcell pressures and saturation-dependent terms are treated as "known", with 
20 such known values being those resulting from the updated linearization. The computed 
saturations can then be used in the next Newton iteration for the timestep. 

Although the additional computations performed in the practice of this invention 
can be applied to all gridcells in the simulation process, preferably the additional 
computations are applied only to those gridcells which incur significant saturation 
25 changes during a timestep. 

BRIEF DESCRIPTION OF THE DRAWINGS 
The present invention and its advantages will be better understood by referring to 
the following detailed description and the attached Figures. 

Fig. 1 (PRIOR ART) is a schematic diagram of the basic steps in a conventional 
30 reservoir simulation process. 

Fig. 2 is a schematic diagram illustrating the basic steps of a reservoir simulation 
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process in accordance with the practice of this invention. 

Fig. 3 is a perspective view of three gridcells illustrating gridcell nomenclature 
used in the description of this invention. 

Fig. 4 is a plot of oil relative permeability (ko) as a function of oil saturation (S 0 ) 
5 in a typical reservoir. 

Fig. 5 is a plot of water relative permeability (k w ) as a function of water 
saturation (S w ) in a typical reservoir. 

Fig. 6 is a plot, similar to Fig. 4, illustrating an example of the linearization of the 
water relative permeability to water saturation relationship. 
10 Figs. 7 A and 7B show a flow chart of various steps used in an embodiment of the 

invention applied with the SEQ method. 

Figs. 8 A and 8B show a flow chart of various steps used in an embodiment of the 
invention applied with the Fully Implicit method. 

1 5 DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

The following description makes use of a large number of mathematical 
symbols, many of which are defined as they occur throughout the text. Additionally, 
for purposes of completeness, a symbols table containing definitions of symbols used 
herein is presented following the detailed description. 

20 Overview 

Before proceeding with the detailed description, a brief description of 
mathematical models and solution methods used in modeling underground 
hydrocarbon-bearing formation will be presented to aid the reader in understanding 
the invention. 

25 In petroleum reservoir simulation, digital computers are used to solve the 

mathematical equations that govern the behavior of fluids in porous media. They 
provide an approximate approach using a suitable gridding format that can 
accommodate any reservoir description. Since persons skilled in the art are familiar 
with such equations and methods used to solve them, it is not necessary to describe 

30 them here. Nevertheless, an overview of general principles of reservoir simulation is 
provided to help those not skilled in the art of reservoir simulation have a better 
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understanding of the invention. 

This overview uses a simplified treatment of the equations. This treatment is 
based on several assumptions that would not normally be valid. First, only two 
substances (oil and water) are assumed to be present, and these fluids are assumed to 
5 form the two phases present. Second, properties of these fluids, such as formation 
volume factor and viscosity, are assumed to be constant. Third, effects of capillary 
pressure are assumed to be negligible. Fourth, the reservoir being modeled is 
assumed to be horizontal. Fifth, flow in the reservoir is assumed to be one- 
dimensional. Sixth, all reservoir properties are assumed to be constant. Seventh, all 
10 cells are assumed to have the same dimensions. 
Material Balance 

In reservoir simulation, the volume of the reservoir is typically discretized into 
thousands of cells or gridblocks. For the purposes of this overview, however, only a 
few cells will be used. Regardless of the number of cells, the numerical model is 
15 based on two fundamental principles. First, each cell must satisfy material balance. 
For oil, this means that the oil in a particular cell at the end of a predetermined 
timestep must satisfy the following equation: 

Oil at New Time = Oil at Old Time + Oil Flow In - Oil Flow Out (1) 

A similar relationship applies for gas, water, and any other fluids in the 
20 reservoir. 

Volume Balance 

The second principle is that each cell must contain the amount of fluid 
required to fill it at a given time. For example, if a reservoir contains only oil and 
water then: 

25 Oil Volume + Water Volume = Cell Volume (2) 

Notation 

Fig. 3 graphically illustrates three cells of a one-dimensional cell system. The 
convention used in this description is to identify cells by means of indices (i, i+1, and 
i-1) as illustrated in Fig. 3. The timesteps (not shown in Fig. 3) are identified by 
30 superscripts as n, n+1, n+2, etc. For example, S " (i is the saturation of the oil phase at 
the beginning of timestep n in cell i, and S**. 1 is the saturation of oil in cell i+1 at 
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the beginning of timestep n+1 (which is the same as the end of timestep n). The 
distance between the centers of two adjacent cells, such as cell i and cell i-fl, is 
represented as Ax. The cross sectional area between the two cells is represented as A. 
Ax, A, and any other quantity that has neither a subscript nor a superscript are 

k 

5 constant. The mobility of oil in cell i is X 0 ± , which is equal to — — , the oil relative 

Ho 

permeability (k rot i) in cell i divided by the oil viscosity (\i 0 ). 
Dairy's Law 

Referring to Fig. 3, flow is assumed to be from left to right; and, as a result, 
cell i is upstream of cell i+1. If upstream mobility weighting is used, as is customary, 
10 mobility at cell i is used in computing flow between the two cells. The rate of oil 
flow from cell i into cell i+1 is represented by the following equation, which is a 
discretization of Darcy's law. 

- ^ikA(pi-p M ) 

F . , = (3a) 



o.i+- Ax 



A similar expression can be written for flow of oil from cell i-1 to cell i. 
15 F , 2iizJ *±L (3b) 



o,i-- Ax 



Relative Permeability Characteristics and Implications 

By definition, mobility (X) is proportional to relative permeability. Fig. 4 is a 
graph of the relative permeability of oil (k r0 ) as a function of oil saturation (S 0 ) in a 
typical oil-bearing formation. The readiness with which oil flows through the 
20 formation decreases as the oil saturation in the formation decreases. The relationship 
between k ro and S 0 depends on the reservoir rock and may vary from one formation to 
another. 

Fig. 5 is a plot of the relative permeability of water (lc™) as a function of water 
saturation (S w ). Similar to the relationship of kr 0 to S 0 in Fig. 4, the readiness with 
25 which water flows through the formation decreases as the water saturation in the 
formation decreases. 

As shown in Fig. 4, oil relative permeability is zero at oil saturations below 
0.2. This means that oil cannot move if its saturation is smaller than 0.2 and, in effect, 
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that the oil saturation cannot become smaller than 0.2. Similarly, Fig. 5 shows water 
relative permeability to be zero at water saturations smaller than 0.25, and water 
saturation cannot become smaller than 0.25. Given that water saturation cannot drop 
below 0.25, oil saturation cannot become higher than 0.75, since S 0 = 1 - S w . Thus, 
5 the valid range of oil saturations is given by 0.2 < S 0 ^ 0.75. Similarly, the valid 
range for water saturation is 0.25 < S w ^ 0.8. The use of upstream weighting helps 
confine saturation to these valid ranges, since a phase can no longer flow out of a cell 
when its relative permeability becomes zero. 
Transport Equations 

10 The amount of oil in a cell is a function of the cell pore volume (V p ), the oil 

saturation (S Q ) in the cell, and the oil formation volume factor (B c ), and can be 
represented by the following formula: 
V S 

Oil in a Cell = -2-^ (4) 

Here B 0 is the oil formation volume factor, which is the ratio of the volume occupied 
15 by oil in the reservoir to its volume at standard surface conditions. 

The pore volume is normally a function of pressure, but it is not being treated 
as such in this overview. In a predetermined timestep, the change in the amount of oil 
in a given cell can be represented as: 

V / x 

Oil at New Time - Oil at Old Time = (S" +1 - S" ) (5) 

20 For a given timestep ( At ), the material balance equation of a cell can be 

represented as 
V V 

— S£ +1 - — S*= At (Oil Flow Rate In - Oil Flow Rate Out) (6) 

This can be rearranged to yield an equation for oil saturation change in a cell 
over the timestep. 

B At 

25 S* +1 - Si - -2 (Oil Flow Rate In - Oil Flow Rate Out) (7) 

V p 

If the preceding equation is applied to cell i, the Oil Flow Rate In is the flow 
rate from cell i-1 to cell i, which is given by equation (4). Similarly, the Oil Flow 
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Rate Out is the rate from cell i to cell i+1, which is given by equation (3). 
Substituting those equations into equation (7) yields 



S 0<i - S o4 - — ( — — - ) (8) 

This equation can be solved for the new oil saturation, resulting in what is 
5 called the oil saturation equation. 

a -._ cn , BpAt A, 0tl .,kAfe> t ., - Pi ) X^kAfa ~ p^,) 

S 0<i = S 0 , i+ ^ ( ^ - ^ ) (9) 

The same equation applies for water, with the subscript o being replaced by w. 
Thus, the water saturation equation is 

Sw ' i = S - i V p ( Ax " Ax ) m 

10 An implication of the volume balance principle, equation (2), is that the sum 

of the water saturation and the oil saturation equals unity, or 

S 0 + S w =l (11) 

If equations (9) and (10) are added and the principle of equation (11) is used, 

the result is 

)kA(p M - Pi ) 

V Ax 

(B 0 X 0 .i + B w X Wii )kA(p, -p,J 
Ax 

This can be simplified to 

+B w X WJ .,)frw -P«)-(B„^ + B w X WJ )(p t -p i+l ) =0 (13) 

This equation can be rearranged to yield the pressure equation used by the IMPES 
20 method and the SEQ method. 
(B o a. 0>i _, +B w X Wii _,)p,_, 

- (B 0 ^o,i-i + B w ^w,i-i + B o X 0 ,i + B w X w> i)Pj 

+ (B o X Oii +B w X Wji )p i+1 = 0 (14) 

Note that no time superscripts appear in this equation. In this overview, the fluids are 
25 assumed to have constant densities. If density was allowed to vary, it would be 
necessary to introduce time into the above equation. Also, recall that, although 



B c At X.i-Mki-i " Pi) KMiPi ~ P i+ i), 
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saturation does not appear in the equation, mobility (X) is a function of saturation. As 
a result, pressure indirectly depends on saturation. 

The SEQ method requires a special form of the saturation equations, which 
can be obtained as follows. The total fluid flow rate is equal to the sum of the oil and 
5 water flow rates. For flow between cells i and i+1, this can be written 



F . , =F . , + F , (15) 

T,H~ o,i+- w,i+- 
2 2 2 

Equation (3) provides an expression for F , , and a similar expression can be 



written for F , . Substituting these into the above equation yields 

w,i+— 
2 

„ ^ 0 ,jkA(Pi-Pi + i) , ^w,ikA(Pi-Pi + i) 

F ,= — ■ + (16) 

T,i+- Ax Ax 

10 This can be solved for the pressure difference, p. - p i+1 , yielding 

Ax 

Substituting this expression into equation (3 a) and canceling certain terms 

yields 

F. ,=-^— F , (18) 

15 This can be written 

F . ,=f 0 ,iF T . , (19) 

where f 0>j is the oil fractional mobility, defined by 



Equation (19) can be substituted into the oil saturation equation (equation (9)) 
20 to yield 

S^SV + ^Af^F , -f 0 .F ' ) (21) 

V p t,- T,, +i 

The equivalent expression for water is 
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S^ l = S:+^(f w ^F , -f wj F ,) (22) 

V p ~2 7 ' + 2 

Solution Methods 

The IMPES method begins each timestep with a pressure solution using 
equation (14). The mobilities are evaluated at the beginning of the timestep, which 
5 means they are based on the saturations at the beginning of the timestep. The 

pressures that are computed are interpreted as applying at the end of the timestep and 
throughout it. The resulting equation is 

( B o^O,i-l + B w^w.i-l )Pm 

10 + (B o X n oJ +B w X n wJ )p? + Y-0 (23) 

Once pressures have been determined, saturations are computed. The 
saturation of either phase can be computed. Assuming that oil saturation is computed, 
the resulting expression is 



„„-, B 0 At Ai-.Wpr-'-pr 1 ) ^(pr'-pff) , 



Sn+1 _ on -^o" * , "0,1-1 yri-i ri / o,i yc-i r iti / . . 

V p Ax Ax 

kAlp" +1 — p" +1 ) 

1 5 Assuming that for whatever reason °' M ^ Fi ; is smaller than 

Ax 

— S ^i^, then S" +1 will be smaller than S" . How much smaller it is will 

Ax 

depend on At. The larger At is, the smaller S^* 1 will be. A value of At can be found 
that will make S" +1 equal to zero, equal to -1, or equal to any other value desired, as 

long as it is smaller than S" . 

20 This is a problem because there is a certain valid range for oil saturation. In 

principle, it can range from zero to one. In actuality, the nature of the relative 
permeabilities restricts it to a narrower range. For the relative permeabilities 
illustrated in Fig. 4., for example, the range is 0.2 < S 0 ^ 0.75. If the timestep ( At ) 
becomes too large, the computed oil saturation can fall below this valid range, a result 

25 that is incorrect. 
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10 



15 



25 



In addition, this behavior can lead to stability problems. On one timestep in 
this example, the saturation drops too low, perhaps below its residual of 0.2. On the 
next timestep, oil cannot move out of the cell because of its low saturation, but it can 
flow in. The saturation becomes too large. Because the saturation is now too large, 
too much oil can flow out, so its saturation drops too much on the next timestep, and 
so on. The remainder of this overview discusses methods intended to address this 
stability problem. 

The SEQ method begins each timestep the same way the IMPES method does, 
by solving a set of equations (23) for pressure. Given these pressures, it next 
computes total flow rates between cells using equation (16) with mobilities evaluated 
at the old time level and pressures evaluated at the new time level. 

Ax Ax 
Given these total flow rates, it is now desired to compute the new saturations 
using equation (21) with the fractional flows evaluated at the new time level. 

Sn+l _ en , 
o,i " a o,i + 



(25) 



? n+l_ on . B e At fn +\ tj 



V ' °- i -' i " T ,i-i- f ".' ,F T, + I ) ' 

P 2 2 



(26) 



However, this cannot be done because the oil fractional flow (£>) terms depend 
on the oil saturations, which are not yet known at the new time level. Instead, these 
terms are linearized with respect to oil saturation, resulting in the following 
expression. 

df" 5 



.(27) 



Substituting equation (27) into equation (26) yields 

j,rn / 
a ^o,i-l 



cn+l c n . B D At 
a o,i - ^o.i "*" w 
P 


B 0 At 


f . + K " 




0> ' dS 0ii 



T,i-- 
2 



.(28) 



T, + - 



With the exception of the new time level saturations, all of the terms in this 
equation are known. A set of equations is formed by writing this one for each cell. 
This set can be solved to yield the new saturations. 
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As described above, the SEQ method is purely sequential. First, pressures are 
computed, then saturations are computed, then the timestep is finished. This can 
cause problems in areas of the formation having high flowrates during a timestep. 

Fig. 6 shows the same graphical relationship of relative water permeability and 
5 water saturation as Fig. 5, except that a line 20 has been drawn tangentially to the k w - 
S w curve 21 at the point on curve 21 that represents a water saturation of 0.4. Line 20 
therefore represents the linearization of the water relative permeability that would be 
used at a water saturation of 0.4. 

Assume that the timestep size is such that 100 pore volumes of fluid pass 

10 through a gridcell in a timestep, the current water saturation in this gridcell at the 
beginning of the timestep is 0.4, and that 10% of the fluid entering the gridcell in the 
timestep is water. This means that the volume of the water entering the gridcell 
during the timestep will be 10 of the gridcell's pore volumes. Clearly the gridcell 
cannot hold this much water, so most of it must flow out during the timestep. 

15 Furthermore, since the current water saturation is 0.4 and the maximum possible 
water saturation is 0.8 (1.0 minus the minimum oil saturation of 0.2), at least 9.6 of 
the 10 pore volumes of water entering the gridcell during the timestep must also leave 
it. In fact, the amount of water leaving the gridcell during the timestep must be 
between 9.6 and 10 pore volumes (assuming that water saturation is increasing). 

20 The water relative permeability must be large enough to permit this to happen. 

Assume that, given the oil relative permeability and oil and water viscosities, a water 
relative permeability of 0.048 would lead to a flow of 9.6 pore volumes and 0.05 
would lead to flow of 10 pore volumes. If the flow into the gridcell is correct, the 
water relative permeability must lie in the narrow range between these two values. 

25 The question is at what water saturation is such a relative permeability 

obtained. The linearization of water relative permeability about a water saturation of 
0.4 is given by 

kJS" =0.008 + 0.12(S W -0.4) (29) 

This expression requires a water saturation of 0.75 to yield a water relative 
30 permeability value of 0.05. As a result, solution of the equations stemming from the 
linearization would yield a water saturation slightly smaller than 0.75, and this value 
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would be used to begin the next iteration. This is significantly in error, since the 
actual function yields 0.05 at a saturation of less than 0.53. It would be much better to 
begin the next iteration with a value around 0.53 than one around 0.75. 

This erroneous result can, in principle, be eliminated by performing the SEQ 
method iteratively in a fashion similar to that of the fully implicit method, which is 
described below. 

The SEQ method is much more stable than the IMPES method, which means 
that it can take much longer timesteps, thereby reducing computing time 
requirements. However, there are some reservoir conditions in which its stability is 
inadequate. The Fully Implicit method is often used in these reservoir conditions. 

The Fully Implicit method simultaneously solves for oil and water saturations 
at the new time level using equations (9) and (10) with all terms evaluated at the new 
time level. 

«■♦■- s- . + B Q At c-,kA(pr-V-p. n+l ) _ o^pr'-pff) . 



V 



Ax 



Ax 



..(30a) 



w.i " S wJ — ( ~ ~ ) 



..(30b) 



However, these are nonlinear equations in which none of the terms can be 
known until all are known. They must be solved iteratively. From the definition of 
the interblock flow, equations (30a) and (30b) can be rewritten 



? n+l_ c n , B pAt 



Sn+1 on , 
o,i " ^o,i + 



7 n+l 
. t 
V °''-2 



-F 1 



n+1 



0,1 + 



•(3D 



2J 



on+l_ c n , B w At 



( 



pn+I _ pn+l 



2 J 



.(32) 



The flow terms are linearized with respect to saturation and pressure, making 
use of the fact that S w = 1 - S 0 . The result is 
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<\n+l,ki_ A /_n+l,k ^n+l,k \ 



F . - 

o,i+^ AX 



^ uk kA(s Pl n+1 ' k -5p^ uk ) 
Ax 



+ ^ ^_ (33) 



+ _2d_ss;* 

dS oJ 0J Ax 
where 

c«n+l,k „n+1,k+l ^n+l,k 

SPi =Pi -Pi (34) 

j?on+l,k c n+l,k+] c n+l,k 

6S o,i = S o,i " S o.i (35) 

and the superscript k denotes iteration number. The equivalent expression for water is 

^n+Uk-i a / n+l,k „n+l,k\ 
F n+l,k+l _ A w,i KAyPi -Pi-H / 

w.i+i Ax 
2 



X»fkA(5pr l *-8pKS l - t ) 



*^ — s < 36 > 



jn,n+l,k UA^ n+1 » k ^n+l.k^ 

dS oJ 0,i Ax 



Substituting equations (33) through (36) into equations (31) and (32) yields a 

set of equations, two at each cell, in which the unknowns are 8S" +1,k and 8I^ n+I,k . 

This set of equations can be solved using methods known to persons skilled in the art. 

10 SS"* 1,k and 5P* 1 +u are solved iteratively until the iteration has converged, at which 

point they will have become negligibly small. Typically, several iterations are 
required. 

A drawback of this method is its per-timestep computational expense. In this 
simplified model, the set of equations that must be solved has two equations for each 

15 cell, rather than just one as does the set of the IMPES method or the SEQ method 
pressure equations. Furthermore, the set of equations must be solved for each 
iteration, of which there are typically several. Overall, a typical computational cost of 
a fully implicit timestep is ten times that of an IMPES method and five times that of a 
SEQ method timestep, but it can be much more. 

20 Furthermore, the method does not always converge. When it fails, the 
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timestep must be repeated with a smaller At. 

Preferred Embodiments 

The inventor has discovered that in solution of the governing flow and 

conservation equations used in simulating fluid properties and behavior in physical 
5 systems, the linearized values of quantities used to compute flow can be more 

accurate than the values of the solution variables from which the linearized values are 

derived. This discovery forms the basis for the simulation process of this invention. 

In the example discussed above with reference to Fig. 6, the linearized value of 

relative permeability differed significantly from the true value. It would normally be 
1 0 expected that the difference constituted an error in the linearization. Surprisingly, in 

reservoir simulation calculations the linearized relative permeability is often more 

accurate than the saturation. 

Although the primary focus of reservoir simulation is the hydrodynamics of 

flow within a reservoir, it also can involve modeling the total petroleum system which 
15 includes the reservoir, the surface facilities that contain fluids, and any interrelated 

significant activity. Those skilled in the art will recognize that the invention can also 

be applied to the equations governing flow through non-reservoir fluid-containing 

media such as wellbore tubing, surface flow lines, separators, and other facilities. 

Fig. 2 illustrates in schematic form the basic steps of reservoir simulation in 
20 accordance with the practice of the present invention. Steps 1 through 4 of Fig. 2 are 
substantially the same as the corresponding steps of a conventional reservoir 
simulation process outlined in Fig. 1, except that in the practice of this invention the 
linearization of the finite difference equations (step 3) uses either the Fully Implicit 
method or the SEQ method. If the Fully Implicit method is being used, step 4 solves 
25 for pressures and saturations simultaneously. If the SEQ method is being used, step 4 
solves first for fluid pressures and then for saturations. The inventor has discovered 
that the results obtained from step 4 can be substantially improved by performing an 
additional computation. After step 4 has been completed, as illustrated in Fig. 2, a 
fifth step of performing a gridcell-by-gridcell computation is made to improve the 
30 prediction of reservoir behavior. First, to prepare for this step, linear estimates of 

nonlinear fluid and transport properties at the end of the time interval are evaluated at 
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the solution (step 4) using the linearizations employed in linearizing the finite 
difference equations (step 3). In the gridcell-by-gridcell computation, the linear 
estimates of fluid and transport properties are used in one or more gridcells adjacent 
to a particular gridcell, while in the particular gridcell, fluid and transport properties 
5 are computed based on the improved solution in that gridcell. 

The step 5 computation can be applied to all gridcells in a grid system being 
modeled. However, to economize on computational time, the additional computation 
is preferably applied only to those gridcells that undergo significant fluid saturation 
changes during a timestep. If the Adaptive Implicit method version of the Fully 
1 0 Implicit method is used as the form of equation solution, it is further preferred that the 
additional computation be applied only within fully implicit regions of the Adaptive 
Implicit method. 

The present invention does not limit the number and type of solution variables 
that can be predicted for a reservoir. Non-limiting examples of three common sets of 
15 solution variables are as follows: 

1 . Pressure, two saturations, selected mole fractions (in models with more than 
two components), and temperature (in thermal models); 

2. Pressure, mass of each component, and energy (in thermal models); and 

3. Pressure, total density of each component, and energy density (in thermal 
20 models). 

Other quantities can be derived from these variables, such as oil, gas, and/or water 
production rates. 
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The first step in the practice of this invention is to equate the formation being 
modeled to a volumetric system comprising a plurality of gridcells. The practice of 
this invention is not limited to any particular type or size of gridcells. For purposes of 
describing the invention only, Fig. 3 illustrates a Cartesian gridcell or gridblock 
5 system. The gridcells may, however, be of any geometric shape, such as 

parallelepipeds (or cubes) or hexahedrons (having four vertical corner edges which 
may vary in length), or tetrahedrons, rhomboids, trapezoids, or triangles. The grid 
may comprise rectangular cells organized in a regular, structured fashion, or it may 
comprise cells having a variety of shapes laid out in an irregular, unstructured fashion. 
1 0 Examples of structured and unstructured gridcell systems are described in Gridding 
Techniques in Reservoir Simulation, by Heinemann, Z. E. and Brand, C. W., lst/2nd 
Stanford University & Leoben Mining University Reservoir Simulation International 
Forum, Alpbach, Austria (Sept 1988/Sept 1989). 

Finite difference equations representative of rock and fluid properties for each 
15 fluid are constructed for each gridcell. The fluid properties in these equations are 
expressed as being implicit in time. Examples of such rock properties include 
porosity, capillary pressure, and relative permeability to each fluid phase. Examples 
of fluid properties include oil viscosity, oil formation volume factor (B 0 ), and initial 
pressure, temperature, and saturation of all fluids in each gridcell. These properties 
20 can be obtained in whole or part from actual reservoir data, or they can be determined 
experimentally or estimated, depending on the type of reservoir simulation being 
performed and the availability of actual reservoir data. Construction of suitable finite 
difference equations representative of reservoir conditions is well known to those 
skilled in the art. 

25 The description of the invention that follows is specific to its application in 

connection with finite difference methods. Those skilled in the art will recognize that 
the invention can also be applied in connection with finite element methods or finite 
volume methods. When it is applied with finite element methods, the gridcells 
become finite elements, and when it is applied with finite volume methods, the 

30 gridcells become finite volumes. The finite difference methods, the finite element 
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methods, and the finite volume methods reduce partial differential equations to a 
finite-dimensional system of algebraic equations. Based on the teachings set forth 
herein with respect to discretization using the finite difference methods, those skilled 
in the art can use the finite element and the finite volume methods in the practice of 
5 the present invention. 

The nonlinear fluid and transport properties in the finite difference equations 
at the end of each time interval are determined by performing an iterative 
computation. For each iteration within this iterative computation, the linear 
approximations of fluid and transport properties are determined, and these linear 

10 approximations are used to derive linear equations. The linear equations can be of the 
form used by either the Fully Implicit method or the SEQ method. The linear 
equations are then solved to yield an estimate of the solution at the end of the time 
interval. In the practice of the present invention, the following additional computation 
is performed to determine an improved solution based on the results obtained in 

15 step 4. 

First, determine at which gridcells the additional computation should be 
performed. As stated above, although the additional computations of this invention 
can be applied to all gridcells, to save computer time, the additional computation is 
preferably applied only to gridcells which are small compared to their volumetric 

20 throughput during a timestep or exhibit large solution changes during a timestep. The 
additional computation can, for example, be applied to gridcells having more than 
about two gridcell pore volumes of water passing through them during a timestep or 
to cells at which saturation changes more than ten percent during a timestep. Persons 
skilled in the art can readily determine which gridcells undergo such changes. As part 

25 of the computations performed in step 4 (Fig. 2), quantitative criteria can be used to 
establish which cells qualify for the additional computation. Non-limiting examples 
of suitable criteria may include pressure change, saturation change, mole fraction 
change, temperature change, residual error at the beginning of the outer iteration, and 
other similar quantities. 

30 Second, compute linear estimates of nonlinear fluid and transport properties 

using the linear approximations of these properties evaluated at the estimated solution. 
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These linear estimates will be needed at all gridcells that are adjacent to, and are. 
therefore in direct-flow communication with, the gridcells at which the present 
invention's additional computation is being performed. If, for example, the additional 
computation is to be performed at gridcell i, the neighboring gridcells may include 
both gridcell i-1 and gridcell i+1, since fluids flow from gridcell i-1 to gridcell i and 
fluids flow from gridcell i to gridcell i+L 

Third, calculate the improved solution at the gridcells in which the additional 
computation is being performed. In this calculation, computation of flow between a 
particular gridcell and the gridcells adjacent to it is based on nonlinear fluid and 
transport properties computed at the improved solution in the particular gridcell and 
the linear estimates of nonlinear fluid and transport properties in the adjacent 
gridcells. This calculation will, almost always, require a nonlinear iteration. This 
iteration can be performed by a variety of methods. All iterative methods involve 
starting with some initial guess for the unknown term and applying some repetitive 
calculation to improve the guess until, after a sufficient number of iterations, the 
simultaneous equations are satisfied to within a predetermined tolerance. The 
solution is said to have converged when calculated values do not change between 
iterations or when the differences fall within an acceptance criterion. It is desirable 
that the iterative method chosen be fast, but it is more important that it be reliable. 

As a last step, the iterative computations are repeated to satisfy a 
predetermined convergence criteria. 
Example 

The following example is provided to further illustrate the practice of this 
invention. In this example, the benefits of performing the additional computation of 
step 5 (Fig. 2) are simulated. It was assumed that oil and water flowed in a one- 
dimensional model similar to the gridcells illustrated in Fig. 3. The oil and water 
relative permeabilities were assumed to have the simple functional forms 



si 



(37) 



and 




(38) 
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Oil and water viscosities and formation volume factors were equal to one. The initial 
oil saturation was constant at 0.9. Water was injected at the left-most gridcell, and the 
injection rate-timestep product was such that ten gridcell pore volumes were injected 
in a timestep. Further details concerning the equations used to simulate this example 
5 are provided in the attached Appendix. Under these assumed conditions, the SEQ 
method yielded the solution set forth in Table 1. 



Table 1 - SEQ Results for Timestep Size of 10 Gridcell Pore Volumes 



Gridcell 
Number 
0) 


1 


2 


3 


4 


5 


6 


7 




-6.2503 


-1.1313 


0.3229 


0.7361 


0.8534 


0.8868 


0.8962 


rUU 


0.7150 


0.9182 


0.9759 


0.9923 


0.9969 


0.9982 


0.9986 . 



Referring to Table 1, S£J/ is the oil saturation in gridcell i computed for the 

10 first timestep (denoted by the first superscript) and the first outer iteration (denoted by 
the second superscript) of this timestep. The superscript £ indicates that this is the 

linearized value. f£]'' is the oil fraction flow defined in a similar fashion. The outer 

iteration will be discussed below. The computed saturations in gridcells 1 and 2 are 
-6.2503 and -1.1313, respectively. These are erroneous by a large amount, since 

15 saturation cannot be less than zero. Table 1 's third row shows the fractional flow 
computed by linearization. These values, though perhaps not correct, are at least 
feasible and are more realistic than the saturations upon which they are based. The 
inventor has discovered that better accuracy can be obtained by determining 
saturations based upon the linearized fractional flows, rather than by the conventional 

20 practice of determining fractional flows based upon saturations. This discovery is 
exploited in two embodiments described below. One embodiment is applied with the 
SEQ method and the other embodiment is applied with the Fully Implicit method. 
The Invention as Applied with the SEQ Method 

The conventional SEQ method makes gridcell-by-gridcell computations to 

25 solve linear equations. At a given gridcell, neighboring gridcell values of the various 
nonlinear functions (such as fractional flow in the example above) are needed. In the 
SEQ method in accordance with this invention, flow into and out of a particular 
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gridcell is deteraiined using linearized values of relative permeability and capillary 
pressure in gridcells neighboring that particular gridcell. Within the gridcell itself, a 
nonlinear, iterative computation is performed. Flow out of the gridcell is determined 
based on the results of this iterative calculation. This iteration is similar to that used 
by the Cascade method. Details of how to perform such iterative computations would 
be familiar to those skilled in the art. 

Consider any gridcell other than the first in the example above. In this simple 
example calculation, the only neighboring gridcell quantity needed to compute its 
saturations is the oil fractional flow at the gridcell immediately left of it. Thus, 
computing saturation at gridcell i consists of solving the problem defined by 



on+l,k+l _ on , B 0 At 



\ 



T-n+l.k+U* _ pn+l,k+l 

r . 1 r . 1 

0,1 — o,i+- 
K 2 2 j 



.(39) 



where 



pn+l,k+l,£ = 



r-n+l.k 
x o,i-l + ' 



df 



n+l,k 



dS, 



o,i-l fcn+l.k+l/ r.n+l,k \ 



o,i-l 



TJ 1 



.(40) 



F n + .,k + l =f n+.,k + l F ( (41) 

o.i +i ' T,, +i 



In equation (39), F n+1, 1 lc+1 

o,i+- 
2 



is a function of f^ y+] 9 which is a function of 



in+Uk+l 
>0J 



As a result, equation (39) is solved iteratively. The particular iterative 



method used is not important, but it is preferably reliable. The iterative method used 
in this example was the reguli falsi technique. Table 2 shows the result of this 
computation. 



20 Table 2 - Saturations Obtained after One Outer iteration when Applying the 
Invention in the SEQ Method 



Gridcell 
Number 
0) 


1 


2 


3 


4 


5 


6 


7 


cU.f 


-6.2503 


-1.1313 


0.3229 


0.7361 


0.8534 


0.8868 


0.8962 




0.1976 


0.5406 


0.6876 


0.7945 


0.8598 


0.8874 


0.8964 
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The computed saturations, , look reasonable. However, they are incorrect. 

During the timestep, 10 gridcell pore volumes of water are injected. As a result, the 
water saturation change totaled over all gridcells must equal 10, and the total oil 
saturation change must equal -10. The initial oil saturation is 0.9. For the 7 gridcells 

5 shown, the total change from this value for the linearized saturations (in the middle 
row of the table) is -9.9862, which is nearly -10. For the computed saturations (in the 
bottom row) the total change is only -1.4361. This is essentially the total change over 
all gridcells, since only very small saturation changes occur in gridcells that are not 
shown. The computed oil saturations indicate that far less oil has been displaced than 

10 water has been injected which is not realistic. 

The newly computed saturations in Table 2 (bottom row) cannot be used as 
they are since they too are erroneous. An outer iteration procedure, similar to the 
Newton method described above for the Fully Implicit method, is applied to the SEQ 
method to obtain more accurate saturation values. Let, for example, S"t u be the oil 

15 saturation after the outer iteration; this is the current estimate of the solution. Let 
gn+i.W b e linearized oil saturation computed during the k th Newton iteration. 
Then the computation is 



Q n+l,k+1 _ on , B Q At 



f \ 

pn+l.k+l.f _pn+l,k+l 
V 2 2 J 



.(42) 



where 



20 F ^U + ^ = f ^^k + M F k + i i (43) 

o.i- ' T.«- 

p n+l,k+l,t _ rn,k , QI o.i-l fon+l,k+l,* Q n+l,k \ /AA , 
1 o,i-l - 1 o,i-l + 17i V^o.i-I ~ ^o.i-l / t 44 ) 

d So,i-i 

pii+l.k+l _ £n+l,k+lpk+I 
2 2 

Each Newton iteration can be performed as follows. 

1 . Compute all nonlinear terms as functions of the current estimate of the 
25 solution. 

2. Compute pressures. 
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3. Compute inter-gridcell total flow rates, F n+, f +1 . 

T,i— 
2 

4. Compute linearized saturations, sj£ ,,k+1,f . 

5. Compute updated linearized values of nonlinear terms, f"f l ' k+1 '' . 

6. Compute new saturations, S" * I,)t+1 , by iteratively solving equation (42). 
5 This iterative procedure is repeated until convergence is obtained. Its first 

four of the six steps are identical to those used in the SEQ method when it is applied 
in iterative fashion. 

Table 3 provides results obtained with the iterative SEQ method described 
above. The procedure converges smoothly and rapidly. After four iterations, the 
10 seven gridcells shown have converged to four significant figures. Given that during 
the timestep ten gridceli pore volumes of water were injected into a system that 
contained very little water, this is very good performance. 



Table 3 - Invention When Applied with Iterative SEQ Method 



Gridceli 
Number 
0) 


1 


2 


3 


4 


5 


6 


7 




-6.2503 


-1.1313 


03229 


0.7361 


0.8534 


0.8868 


0.8962 


oU 


0.1976 


0.5406 


0.6876 


0.7945 


0.8598 


0.8874 


0.8964 






0.1976 


0.2204 


-0.2728 


-1.3946 


-1.8423 


-0.6745 


0.3259 


sS 


0.1976 


0.2542 


0.2949 


0.3509 


0.4426 


0.5637 


0.6857 






0.1976 


0.2542 


0.2930 


0.3249 


0.3522 


0.3147 


0.0179 


oU 


0.1976 


0.2542 


0.2930 


0.3238 


0.3498 


0.3729 


0.3964 




a 1.4,* 


0.1976 


0.2542 


0.2930 


0.3238 


0.3499 


0.3730 


0.3942 


nl.4 

^0,1 


0.1976 


0.2542 


0.2930 


0.3238 


0.3499 


0.3730 


0.3942 



1 5 The Invention as Applied with the Fully Implicit Method 

The practice of the present invention can also be applied in connection with 
the Fully Implicit method. Equation (42) is replaced by the following modification of 
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equation (28). 

qn+ u +1 «. , B a At c:r u kAk-Y' k+l -pr u+i ) 

h>_J^ki — iPi±^J) (46) 

Ax 

where 

jn,n,k 

c ^n+l.k+I/ ^n.k ( qA o,i-l f c n+l,k+l,f on+l,k J 

5 *-o,i-l = ^o,i-l + ^ l b o,i-l - ^o.i-l ) (47) 

Each Newton iteration is performed as follows. 

1 . Compute all nonlinear terms as functions of the current estimate of the 
solution. 

2. Compute pressures ( p f +1,k+1 ) and saturations ( S"+ ! ' k+I ' ) . The procedure 
10 used is that for a conventional fully implicit Newton iteration. 

3. Compute updated linearized values of nonlinear terms, A^ Uk+u . 

4. Compute new saturations, S"} , k+1 , by iteratively solving equation (46). 
This iterative procedure is repeated until convergence is obtained. Its first two 

steps are identical to those used in the Fully Implicit method. 

15 Persons skilled in the art will readily understand that the present invention is 

computationally intense. Accordingly, use of a computer, preferably a digital 
computer, to practice the invention is virtually a necessity. Computer software for 
various portions of the simulation process are commercially available (for example, 
software is commercially available to develop gridcells, display results, calculate fluid 

20 properties, and solve linear sets of equations). Computer software for other portions 
of the invention could be developed by persons skilled in the art based on the 
teachings set forth herein. Flow diagrams of a computer program are generally 
illustrated as Figs. 7 and 8. Fig. 7 illustrates a generalized flow diagram of the present 
invention when applied to the SEQ method and Fig. 8 illustrates a generalized flow 

25 diagram of the present invention when applied to the Fully Implicit method. 
Generalization of the Computation 

Define X as the vector of solution variables. In the example used above, 
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So 
Sw 



(48) 



10 



15 



In other models, X can include component masses, mole fractions, 
temperature, energy, and other quantities. In any case, X comprises a set of state 
variables that is sufficient to define the state of the system. Furthermore, it is the set 
of solution variables, and other variables are considered to be functions of X, either 
directly or indirectly. 

The flow between gridcells i and i+1 can be written as 



This equation is written for flow between gridcells i and i+1. An equivalent 
expression can be written for flow between gridcells i and i-1, or in a 
multidimensional model for flow between gridcell i and any other gridcell with which 
it is connected. 

The linearization of the flow between gridcells i and i+1 can be expressed as 



The difference between the iterative SEQ method and the Fully Implicit 
method is in the definition of F. For either method, F includes the effects of pressure, 
saturation, capillary pressure, composition, and anything else that affects flow 
between gridcells and is included in the flow model. 

Two new sets of variables can be defined. Let the vector Y consist of the 
components of X that are used directly in computing flow. The only such component 
may be pressure. Let the vector Z consist of quantities that are functions of X and are 
used in computing flow. Relative permeability and capillary pressure are two such 
quantities. 




(49) 
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Using these, we can rewrite equation (49) as 



F ,=G .(Yj.Z^Y^.ZjJ (51) 

v,i +i v, +i 

and equation (50) as 



5G n + l.k 



.1 -1 2V * ' ' ' 



aG n + l,k 

v,i+- 



dZ: K 1 1 } 

(52) 

5G n+, ' k 



V> 2 Km+l,k+1,£ v n+l,k| 

k 
\_ 

2 (^n+\ t k+\ t e 7 n+l,k \ 
— V^i+l / 



aG n + l,k 

v,i+- 



5 where the linearized value of Zi is defined as 

^n+u+i.r _ z »+u i fj<i+u+u _ Jfc r . n+I '* J (53) 

It would be ideal if the linearization of equation (51) could then be written 

F n + I.k + U = q ( (Y n + l.k + l.^ Z n + l.k + l.^ Y ^ Ml \Z^ MU ) (54) 

Then the invention would be able to use the expression 
10 F n+I * +1 =G , (y* + *+\Zp w 9 Y$*+ u 9 Z$*+ l ' i ) (55) 

in the case where the invention is being applied to gridcell i. Unfortunately, equation 
(55) contains cross-product terms (i.e., terms containing a product such as 
( s n+u + i _ gn+u )^n+i.k+i _ p «+u )) that are missing from equation (52). This 

introduces an consistency problem if equation (55) is applied directly. It is important 
15 that consistency be maintained, which means that steps must be taken to ensure that 
pn+ijc+i =F n + i.k + u when Y ; +1 ' k+I = Y" +I,k+I,f and Z^ 41 = Z? +u+1 '' . This is 

v,i+- v,i+- 
2 2 

accomplished by replacing equation (55) by 
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Fn+Kk+1 (•\rn+hk+\ ^n+Kk+l \rn+Mc+l/ ^n+l.lc+l.f \ rn+l.k+l/ 

i =(j . il Y i »A > Y i + i r h . \ ( 56 ) 

V,l+- V.1+- V.1+- 

2 2 2 

where e b+1 * +1 '' is an error term given by 



i 

2 



gn+i,k+i,*_Q ^Y n+1 ' k+I ^ ^p+i.k+i,* Y n | 1,k+1,f Z"Y* k+1,? )~ F n+lfk+1 *' (57) 

v,i+- v,i+- v,i+- 

2 2 2 

and F n+1, , k+1/ is computed using equation (52). The error term can be evaluated 

v,i+- 
2 

5 based on the linearized values, all of which are known when the invention is applied. 
During application of the invention to gridcell i, equation (56) is used 
regardless of whether gridcell i or i+1 is upstream. If i+1 is upstream for all phases 
and there is no capillary pressure, flow between i and i+1 depends only on the 
linearized values at i+1. If i is upstream, the same flow depends only on the new 

10 values being computed at i. In either case, the computation is the same as it was in 
the examples discussed above. If there is countercurrent flow, equation (56) 
appropriately uses a mixture of linearized i+1 values and new i values. Similarly, if 
capillary pressure effects are introduced, equation (56) uses the appropriate values 
from the two gridcells. 

15 It is expected that the process of the present invention will play an important 

role in lowering the cost of production of oil and gas resources. The economical and 
efficient production of reserves is crucial to prolonging the world's oil and gas 
resources. It is believed that the present invention can make a significant contribution 
to that effort. 

20 The invention is not to be unduly limited to the foregoing which has been set 

forth for illustrative purposes. On the contrary, a wide variety of modifications and 
alternative embodiments will be apparent to persons skilled in the art without 
departing from the true scope to the invention, as defined in the claims set forth 
below. 

25 

Symbols Table 

The most commonly used symbols in this description of the invention are 
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defmed as follows: 

A - cross sectional area between two gridcells 

B 0 = oil formation volume factor (reservoir volume of the oil phase per 
standard volume of oil) 
5 B w - water formation volume factor 

E = error in computed flow between gridcells defined by equation (55) 
f G = fraction of total mobility that is oil mobility 
f w = fraction of total mobility that is water mobility 
F D = rate of flow of oil 
10 F w = rate of flow of water 

Ft = rate of flow of oil plus water 

G = generalized functional form of flow between gridcells 

i = subscript denoting gridcell 

k = permeability; as a superscript denotes iteration 
15 k ro = relative permeability of oil 

k rw = relative permeability of water 

£ - superscript indicating linearized value 

m = superscript denoting inner iteration number 

n = superscript denoting time level; n=0,l,2,3 ... 
20 p = pressure 

q 0 = oil flow rate 

q w - water flow rate 

= oil saturation at current time 
S n = water saturation at current time 

w 

25 Sq +1 = oil saturation at new time 

S n+1 = water saturation at new time 

w 

T = subscript indicating total 

T 0 = gridcell pore volume divided into the sum of oil in a gridcell at the 
beginning of a timestep and the linearized amount of oil flowing in 
30 during the timestep 
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At = change of time 

(i o - oil viscosity 

H w = water viscosity 
V p = gridcell pore volume 
5 Ax = distance between centers of two gridcells 

x - distance 

X = vector of solution variables; see equation (48) 

Y = vector containing the components of X that are used directly in computing 

flow between gridcells 
10 Z = vector containing quantities that are functions of X and are used in 

computing flow between gridcells 
X o = mobility of oil (oil relative permeability divided by oil viscosity) 
X w = mobility of water (water relative permeability divided by water 

viscosity) 

15 8 = operator denoting change over an iteration 

A = operator denoting change in space or time 
v = subscript denoting an arbitrary phase 
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APPENDIX 

SEQ METHOD EXAMPLE COMPUTATIONS 

The invention was tested on the simple one-dimensional series of rectangular 
gridcells in which water is injected into the left block. Figure 3 shows three gridcells 
5 in the series of gridcells. Water, oil, and total mobilities are given by 

X w =Si=(l-Sj (A-l) 

^o=S^ (A-2) 

K =^«+K (a-3) 

Oil saturation is solved for. Water and oil mobility derivatives with respect to 
10 oil saturation are 
dX„ 



dS 0 



= -3(l-S 0 ) 2 (A-5) 



(A-6) 



The oil fractional flow and its derivative with respect to oil saturation are 



f„=^- 



(A-7) 



15 



dS„ 



dX 0 dX„ 

A.,„ — — A. 



dS„ 



dS„ 



The linearized flow during a timestep from gridcell i to i+1 is given by 



* 1 "^v 
oi+- 
2 



jjrn+I.k 
rn+l,k , ai oi con+l»M 



(A-8) 



(A-9) 



where q w is the water injection rate and 8S** l *' { is defined in equation (35). 

Only water, and thus no oil, flows into the first gridcell. As a result, the linear 

20 saturation change in it is 



6S oi = s oi- s oi +- rr~ 



j f n+l,k 



dS 



n+l,k,* 



ol 



.(A-10) 
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Solving this yields 



Xon+l^ - 12 

5S oi " "TTTTnTTk (A-ll) 

t j Qw At df o1 

V P dS 0l 

The linearized fractional flow out of the first gridcell is then computed as 

n+l,k 

f oi = f oi — 5 S oi (A-12) 

where i = 1 . 

Then, given the linearized fractional flow at an arbitrary gridcell i-l, the 
linearized saturation change at gridcell i can be computed by 

r*n c n+1,k , q w At [pn+l.k+U fn+l,kl 
^oi ^oi + ~T7 L r oi-l ~ I oi J 

q w At df oi 



1 + - 



Vp dS„; 



' P w o\ 

The procedure for computing the linear saturation changes and linearized 
10 fractional flows is as follows. 

1 . Compute the linear saturation change at gridcell 1 using Eq. (A-l 1). 

Update the "linear" saturation by S^ +u+1 ^ = S£ Uk + dSlf M . 

2. Compute the linearized fractional flow from gridcell 1 to gridcell 2 using 
Eq. (A-12) with i = L 

15 3. Beginning at gridcell 2, for each gridcell first compute the linear saturation 

change using Eq. (A- 13) then compute the linearized fractional flow using 
Eq. (A-12). Continue until these computations have been performed for 
all gridcells. 

This simplified procedure can be used because the model is one-dimensional, 
20 incompressible, and all flow is in the direction of increasing i. Because the model is 
one-dimensional and incompressible, the total fluid velocity is constant and it is not 
necessary to solve for pressure. Because the model is one-dimensional and flow is 
uni-directional, the computational procedure can proceed from small i to large i, as 
described above. In a more general model, it is necessary to solve matrix equations to 
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first determine pressure changes and then determine linear saturation changes. 

The next step is to compute the improved saturation estimates. The 
computation is performed a gridcell at a time and is based on the following equation. 

c n+l,k+l c n , Qw^t f f n+l,k+U x*n+l,k+l /A , A . 

S oi = S oi + ~r7~L t oi-l ~*oi ( A * 14 ) 

5 Flow into the gridcell is based on the linearized fractional flow. Flow out is 

based on the gridcell's final saturation. The computation must be performed 
iteratively. The reguli falsi method was used. The method works by bounding the 
solution, then repeatedly narrowing the bounds. The following objective function is 
used for the m tb iteration of the reguli falsi method*. 

1A TT m _on+l,k+l,m , Qw^t f n+l.k+l,m T n+Uc+l,£ , A , <A 

AU r obj,i - ^o\ *~r; I oi ~~ A oi 

V P 

where 

T oi = S °l + ~^~~ (A-16a) 
V P 

r n+l,k+l^_Qn ,Qw^ f n+l,k+l^ 



oi = S oi f oi-l > l> 1 (A-16b) 

is the gridcell pore volume divided into the sum of oil in the gridcell at the beginning 
15 of the timestep plus the (linearized) amount flowing in during the timestep. 

The objective of the iteration is to find a value of S£ lMhm that causes F™ ; to 

become zero. The initial bounds on the solution are 

Sn^O • (A-17) 

and 

20 =1 (A-18) 

At these saturations, the objective function is respectively 

F^—Tr^'' (A-i9) 

and 

F max- 1 + -T^ hi .-(A-20) 

V P 

25 If the objective function were linear in saturation, it would be zero at the 

following saturation. 
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Sn+l.k+l,m+l _ o t? ^max ^min 
oi — ^min ~ r min ^ ~ (A-21) 

"max **min 

This saturation is used as the next estimate. The fractional flow is computed at this 
saturation, then the updated objective function, F^*/ , is computed using Eq. (A- 15). 
If the objective function's absolute value is small enough (in the examples in the text, 
5 less than 0.0001) the iteration is terminated. If not, the bounds are narrowed by the 
following procedure. 

IfF<S <0,set 

S min =S- 1MUm+1 (A-22) 

(A-23) 

10 IfF™ 1 >0,set 

Sn.x^S^^^ (A-24) 

(A-25) 

The process beginning with Eq. (A-21) is then repeated until Jf^/J < 0.0001 . 

The outer iteration's final saturation, S"* u+l , is set to the final s^' k+Km+1 . This 
1 5 saturation is used to begin the next outer iteration. 

A timestep is performed in the following major steps. 

1. Initialize the computation by setting S"; equal to sjj* uk+1 , the final 
linearized saturation from the preceding timestep. 

2. Compute the linear saturation changes, saturations, and fractional flows 
20 using the procedure described following Eq. (A-13). 

3. Test for convergence based on the linear saturation change. If the largest, 
in magnitude, linear saturation change is smaller than a pre-specified 
tolerance, convergence has been obtained. The solution for the timestep is 
the set of linear saturations from step 2. 

25 4. If convergence has not been obtained, compute improved saturations using 

the process described beginning with Eq. (A- 14). 
5. Repeat steps 2 through 4 until convergence is obtained. 
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Table 1 presents linear saturations and fractional flows for the first outer iteration of 

the first timestep of a problem with = 10 and S° = 0,9 . Table 2 presents the 

same linear saturations plus the improved saturations obtained in step 2 immediately 
above, again for the first iteration of the first timestep. Table 3 presents the same 
5 quantities, but for the first four outer iterations of the first timestep. 
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WHAT IS CLAIMED IS: 

1 . A process for predicting a characteristic of a physical system containing 
multiple fluids, comprising the steps of: 

(a) equating the physical system to a volumetric grid system comprising a plurality 
5 of gridcells; 

(b) for each gridcell, constructing nonlinear governing equations representative of 
the behavior of the fluids in the gridcell over a time interval, said equations using 
fluid and transport properties computed at the end of the time interval; 

(c) determining linear expressions for nonlinear fluid and transport properties in the 
10 governing equations and using these linear expressions to derive linear equations 

over the time interval; 

(d) solving the linear equations to determine an estimated solution at the end of the 
time interval; 

(e) determining linear estimates of the nonlinear fluid and transport properties at the 
1 5 end of the time interval by evaluating the linear expressions for the fluid and 

transport properties at the estimated solution; 

(f) determining an improved solution by performing iterative gridcell-by-gridcell 
computations wherein fluid flows to and from a particular gridcell are 
determined using (i) fluid and transport properties for that particular gridcell 

20 evaluated at the improved solution in that particular gridcell and (ii) the linear 

estimates of fluid and transport properties at one or more gridcells adjacent to the 
particular gridcell; 

(g) using the results of step (f) to predict a characteristic of the physical system and 
the fluids it contains at the end of the time interval; and 

25 (h) performing steps (b) through (g) for a plurality of time intervals and using the 

results to predict a property of the physical system and the fluids it contains as a 
function of time. 
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The process of claim 1 wherein the physical system is a hydrocarbon-bearing 
subterranean formation. 

The process of claim 1 wherein the physical system comprises fluid- 
containing facilities associated with production of hydrocarbons from a 
subterranean hydrocarbon-bearing formation. 

The process of claim 1 wherein the results of step (f) are used to predict the 
saturation of a fluid in the physical system. 

The process of claim 1 further comprising after step (f) and before step (g) 
repeating the iterative computations of steps (c) through (f) to satisfy a 
predetermined convergence criterion. 

The process of claim 1 wherein the gridcells are finite difference gridcells and 
the governing equations are finite difference equations. 

The process of claim 1 wherein the gridcells are finite elements and the 
governing equations are finite element equations. 

The process of claim 1 wherein the gridcells are finite volumes and the 
governing equations are finite volume equations. 

The process of claim 1 wherein the step (f) is applied only to gridcells in 
which significant changes in fluid properties occur during a timestep. 

The process for predicting behavior of a subterranean formation containing 
multiple fluids, comprising the steps of: 

(a) equating said formation to a volumetric grid system comprising a 
plurality of gridcells; 

(b) for each gridcell, constructing nonlinear finite difference equations 
representative of rock and fluid properties in each gridcell for each 
fluid, said fluid properties being represented as implicit in time; 

(c) linearizing the nonlinear terms in the finite difference equations to 
derive the linear equations used in the sequential implicit method 
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which will include linear pressure equations and linear saturation 
equations; 

(d) solving the linear pressure equations to determine an estimated 
pressure solution; 
5 (e) computing total flow rates between gridcells; 

(f) using the results of step (e) to modify the linear saturation equations; 

(g) solving the linear saturation equations to determine an estimated 
saturation solution; 

(h) improving the estimated saturation solution by performing iterative 
10 gridcell-by-gridcell computations wherein flows to and from a certain 

gridcell are determined by linearized values of nonlinear terms in 
gridcells adjacent the certain gridcell and new timestep values of these 
nonlinear terms for the certain gridcell; 

(i) computing phase flow rates between gridcells based on the results of 
15 step(g); 

(j) computing the amount of each fluid in each gridcell using the results of 
step (i); 

(k) repeating the iterative computations of steps (c) through (j) to satisfy 

predetermined convergence criteria; 
20 (1) using the results of step (k) to predict a property of the formation and 

the fluids it contains at the current time; and 
(m) performing steps (a) through (1) for a plurality of timesteps and using 

the results to predict a property of the formation and the fluids it 

contains as a function of time. 
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